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“MULTISCALE  DYNAMICS  AND  INFORMATION  IN  DATA  COLLECTION  AND 
ASSIMILATION  FOR  ENVIRONMENTAL  APPLICATIONS” 


PI:  N.  SRI  NAMACHCHIVAYA 


This  is  the  final  report  for  AFOSR  Grant  FA9550-12-1-0390,  which  was  awarded  in  07/01/2012.  This 
research  involves  the  work  of  the  PI  and  graduate  students  at  University  of  Illinois  at  Urbana-Champaign 
(UIUC)  in  collaboration  with  a  graduate  student  of  Humbolt  University  (HU),  Berlin. 

1.  Summary  of  Project  Activities  and  Findings 

Data  assimilation  or  filtering  involves  blending  information  from  observations  of  the  actual  system 
states  with  information  from  dynamical  models  to  estimate  the  current  system  states  or  certain  model 
parameters.  The  filtering  problem  relies  on  three  fundamental  ingredients,  namely  1)  sensor  placement: 
where  the  sensors  are  placed  in  order  to  obtain  the  most  useful  information,  2)  sensor  fusion:  how  to 
combine  the  measurements  from  different  sensors,  and  3)  estimation:  how  to  use  the  measurements 
to  obtain  the  best  possible  state  estimates.  In  this  project,  we  considered  the  data  assimilation  prob¬ 
lem  for  multi-timescale  systems.  An  understanding  of  how  scales  interact  with  information  led  to  the 
development  of  rigorous  reduced-order  data  assimilation  techniques  for  these  high-dimensional  problems. 

The  main  results  of  the  research  project  are: 

(1)  Rigorous  mathematical  development  of  a  reduced-order  particle  filtering  method  for  high-dimensional, 
multiscale  random  dynamical  systems. 

(2)  Development  of  a  particle  filtering  method  adapted  to  high-dimensional,  multiscale,  chaotic 
systems 

The  research  has  led  to  three  journal  publications  [1,  2,  3],  four  conference  papers  [4,  6,  7,  11],  and 
two  key  note  lectures  one  at  an  international  research  workshop  [5]  and  the  other  at  an  international 
conference  [10].  Two  Ph.D  students  who  were  involved  in  the  research  are  in  the  course  of  their  graduate 
studies  at  UIUC,  and  one  graduate  student  who  collaborated  in  this  project  has  completed  his  studies 
at  HU,  Berlin. 

1.1.  Reduced-order  filtering.  Dynamical  models  for  describing  climate  evolution  consist  of  coupled 
ocean  and  atmospheric  models,  which  are  high-dimensional  with  multi-timescale  dynamics,  described  by 
stochastic  differential  equations  of  the  form: 

dXet  =  b(XtE,  Zet)dt  +  a{Xet,  ZI)dWt ,  X^  =  xG  Rm, 

dZ$  =  4  f{X$,  Zf)dt  +  -g{Xl  Zf)dVt,  Z%  =  z  £  Rm, 

£-  £ 

where  W  and  V  are  independent  Wiener  processes  that  model  the  uncertainties  in  the  dynamical  system 
and  e  <<  1  is  a  small  parameter  that  characterizes  the  timescale  separation  between  the  fast  and  slow 
processes.  In  (1),  the  fast  and  slow  processes  are  Z  and  X ,  respectively. 

In  the  filtering  framework,  the  objective  is  to  obtain  the  best  estimate  of  (A6,  Ze)  from  observations 
Y,  which  are  described  by 

YtE  =  [T  h(XI,  Zf)ds  +  Bt,  Yq  =  0. 

JO 

Formally,  the  objective  is  to  find  a  normalized  filter  7if : 

tt/M  =  E  y$) ,  y*  =  a  (Yf  :  0  <  a  <  t)  V  V, 

where  <p  €  Rm  x  R™  — >  R  is  any  bounded,  measurable  function  and  yt  is  the  algebra  generated  by  the 
observation  Y  up  to  time  t. 

Data  assimilation  for  high-dimensional  dynamical  systems  such  as  that  described  in  (1)  are  faced 
with  the  computational  complexities  that  arise  from  dimensionality  issues.  Part  of  the  work  performed 
under  this  research  grant  has  been  the  development  of  a  reduced-order  particle  filtering  technique  that 
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enables  data  assimilation  on  a  high-dimensional  system  by  making  use  of  an  associated  reduced-order, 
stochastically-averaged  system. 

In  studying  multiscale  systems  such  as  (1),  we  are  usually  interested  in  the  statistics  of  the  long¬ 
term  (“slow”)  dynamics.  Hence,  we  can  take  advantage  of  the  timescale  separation  between  the  fast 
and  slow  processes  and  the  appropriate  stochastic  properties  that  ensure  the  existence  of  a  homogenized 
system  associated  with  the  original  multiscale  system  (1).  Under  appropriate  conditions  on  the  coefficient 
functions  of  (1),  it  is  known  that,  for  generator  £E  associated  with  (1)  and  probability  density  pe , 

£  0  =>  (ps,  £e)  — >  (p° ,  £), 

where  p°  is  a  reduced  dimension  density  function.  We  can  then  find  a  process  associated  with  £  governed 
by  SDE  of  the  form: 

(2)  dX°  =  b{Xg)dt  +  a(X°)dWt,  Xg  =  x  €  Mm, 


such  that  the  process  Xs  converges  weakly  to  X°  in  the  limit  as  e  — >  0.  The  fast  process  dynamics  is 
averaged  into  the  homogenized  system  (2),  thus  the  dynamics  of  X°  effectively  captures  the  distribution 
of  the  coarse-grained  dynamics  of  (1).  The  idea  behind  the  development  of  a  reduced-order  filter  for  the 
multiscale  systems  is  to  make  use  of  the  existence  of  the  process  (2).  Formally,  the  task  is  to  show  that 
7 re  is  close  to  n°  in  Lp-sense  as  e  — >  0,  where  7r°  is  the  filter  associated  with  the  homogenized  process 
(2).  For  purposes  of  filtering  applications,  this  means  that,  for  systems  with  large  timescale  separation, 
a  reduced-order  filter  ir°  can  be  used  in  place  of  the  original  7re  to  estimate  the  coarse-grained  dynamics. 
This  dimensional  reduction  of  the  filtering  problem  translates  to  reduction  in  computational  costs. 

The  important  result  of  the  analysis  part  of  the  project  in  the  development  of  the  reduced-order 
filtering  technique  is: 


(3) 


lim  supE  [d  (7tf ,  7r°) ]  =  0,  VT  >  0 

t<T 


The  analysis  and  result  for  a  one-dimensional  system  with  timescale  separation  is  shown  in  [12].  In  the 
one-dimensional  case,  a  stochastic  partial  differential  equation  (SPDE)  driven  by  real  observation  was 
constructed  for  the  unnormalized  conditional  density  of  the  reduced-order  filter.  Convergence  was  shown 
by  showing  that  the  error  between  the  solution  of  the  Zakai  equation,  which  governs  the  unnormalized 
conditional  density  of  the  multiscale  system,  and  solution  of  the  SPDE  corresponding  to  the  reduced-order 
filter  vanishes  as  £  — >  0.  The  calculation  required  an  application  of  Levy’s  Theorem  for  time-changed 
representation  of  Brownian  motion,  which  is  holds  only  in  one  dimension. 

In  the  research  under  this  grant,  work  was  done  on  extending  the  result  to  the  general  multi¬ 
dimensional  case.  A  SPDE  driven  by  real  observation  was  constructed  for  the  unnormalized  reduced-order 
filter.  Solution  to  this  SPDE  is  a  measure- valued  process,  and  the  task  was  to  show  that  the  error  relative 
to  the  solution  of  the  measure- valued  Zakai  equation  for  the  multiscale  system  is  bounded  and  vanishes 
as  e  — >  0.  This  was  achieved  by  invoking  dual  representations  of  the  unnormalized  multiscale  filter,  vs , 
and  the  reduced-order  filter  driven  by  real  observation,  v°.  The  dual  processes  satisfy  backward  stochas¬ 
tic  partial  differential  equations  (BSPDEs)  that  are  related  to  the  forward  SPDEs  of  the  corresponding 
filters.  The  outline  of  the  convergence  proof  is  as  follows: 

ve  is  asymptotically  expanded  into  an  order  1  term,  a  correction  term  of  order  e,  and  a  remainder 
order  higher  than  e.  These  expansion  terms  are  governed  by  BSPDEs,  and  the  BSPDE  of  the  order 
1  term  corresponds  to  that  of  v°.  Hence,  the  correction  and  remainder  terms  correspond  to  the  error 
between  ve  and  v°.  Using  BSPDE  theory,  probabilistic  representations  of  the  corrector  and  remainder 
are  obtained  in  terms  of  solutions  of  backward  “doubly”  stochastic  differential  equations.  From  there, 
bounds  of  order  yfe  were  obtained  for  the  corrector  and  remainder  using  existing  probabilistic  estimate 
results  in  literature.  Hence,  the  error  between  ve  and  v°  vanishes  as  e  — >  0.  This  convergence  was 
translated  back  to  convergence  of  the  unnormalized  filters  by  the  dual  relations,  and  boundedness  of  the 
measure  change  (in  the  Kallianpur-Striebel  formula  relating  normalized  and  unnormalized  filters)  gives 
convergence  of  the  normalized  filters.  This  proof  was  completed  at  the  end  of  the  first  year  research. 
A  preliminary  result  was  presented  in  a  conference  and  published  in  the  proceedings  [4],  and  the  final 
results  have  been  published  in  [1]. 


1.2.  Reduced-order  particle  filter  adapted  to  chaotic  systems.  Based  on  the  analytical  re¬ 
sults  previously  described,  a  reduced-order  filtering  algorithm,  the  Homogenized  Hybrid  Particle  Filter 
(HHPF)  was  developed,  described  in  detail  in  [13]  for  the  continuous  time  case. 
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In  the  second  year  under  the  research  grant,  work  was  done  to  apply  the  HHPF  algorithm  to  high¬ 
dimensional  chaotic  systems.  Due  to  the  sensitivity  to  initial  conditions  nature  of  chaotic  systems,  the 
regular  HHPF  was  found  to  be  inadequate  for  estimating  the  system  states  accurately  and  efficiently. 
The  discrete-time  filtering  problem  was  considered  in  the  attempt  to  adapt  the  reduced-order  filtering 
algorithm  to  chaotic  systems.  A  discrete-time  version  of  the  HHPF  was  developed,  which  incorporates 
the  technique  of  Sequential  Importance  Sampling  (SIS)  for  particle  filtering.  The  idea  behind  importance 
sampling  is  to  represent  a  target  probability  density  that  is  difficult  to  sample  from  using  an  alternate 
proposal  density.  The  discrepancy  from  not  using  the  true  density  is  accounted  for  by  appropriately 
weighting  the  proposal  density. 

In  the  filtering  problem,  the  objective  of  importance  sampling  is  to  represent  the  posterior  density  of 
the  state  given  observations  using  a  proposal  density.  Using  SIS,  a  proposal  density  can  be  generated 
using  appropriately  weighted  particles,  where  particle  weights  are  updated  sequentially.  For  the  discrete¬ 
time  case,  the  optimal  proposal  density  (that  minimizes  the  variance  of  particle  weights)  at  a  discrete 
timestep  k.  qopt(xk\xk-i,yk),  is  a  known  result  that  is  the  normalized  product  of  the  prior  density  based 
on  signal  dynamics  and  the  likelihood  of  observation  given  the  state,  with  corresponding  sequential 
particle  weight  update  expression.  Consider  the  discrete-time  signal  dynamics  and  observation 

(4)  Xk  =  B(Xfc_!)  +  axAWk, 

Yk  =  H{Xk)  +  ayAVk. 

For  the  case  of  linear  observation  function  H , 

(5)  qopt(xk\xk-i,yk)  =Af(B(xk-i)  +  a(xk-i,yk),Q), 

Q=  (Q-1  +  HTR~1H)~1 , 
a(xk-i ,yk)  =  QHT R~1(yk  -  HF(xk- 1)), 

where  Q  =f  crxa^ ,  R  =f  UyCTy  .  The  proposal  density  (5)  can  be  generated  by  propagating  particles 
using  an  augmented  signal  dynamics.  For  particle  i.  the  original  dynamics  (4)  is  modified  into 

(6)  =  B{Xk_x)  +  a(Xk_i,  Yk)  +  axWk 

where  ax  is  such  that  <Jxax  =  Q.  Then  Xk  behaves  like  a  particle  sampled  from  gopt(-|x[._1,  yk). 

A  similar  augmented  signal  dynamics  for  particle  propagation  result  was  also  obtained  by  considering 
the  discrete-time  problem  using  a  stochastic  optimal  control  approach.  The  problem  setup  was  to  consider 
particle  dynamics  governed  by  the  stochastic  difference  equation  of  the  form  (6),  where  a(X'k_1,  Yk)  is 
now  an  unknown  “control”  that  is  to  be  applied  to  steer  particles  towards  the  true  signal  based  on 
observations.  This  is  a  one-step-ahead  control  problem.  The  control  term  was  obtained  by  minimizing 
a  quadratic  cost  function  that  penalizes  the  amount  of  control  applied  over  a  timestep  and  the  distance 
of  a  particle  from  the  location  prescribed  by  the  observations.  This  augmentation  of  signal  dynamics  for 
particles  propagation  can  be  related  to  existing  techniques  widely  used  in  geophysical  data  assimilation 
such  as  the  4D-VAR  technique. 

By  implementing  SIS  with  the  optimal  proposal  density  in  the  HHPF,  we  have  a  discrete-time  version 
of  the  HHPF  that  is  adapted  to  systems  with  high  sensitivity  to  initial  conditions.  For  numerical 
experiments,  the  optimized  HHPF  was  applied  to  a  396  dimensions  (36  slow,  360  fast)  Lorenz  ’96  model 
that  mimics  midlatitude  dynamics  of  an  atmospheric  variable.  The  Lorenz  ’96  model  is  a  popular  testbed 
for  numerical  statistical  analysis.  The  purpose  of  the  numerical  experiments  is  to  study  the  performance 
of  the  optimized  HHPF  in  estimating  the  slow  state  variables  based  on  complete  and  incomplete  (not 
every  dimension  of  the  true  signal  is  observed)  noisy  observations  of  the  true  signal.  The  optimized  HHPF 
is  compared  with  the  regular  HHPF  (uses  prior  as  proposal),  an  existing  weights  variance-reduction 
particle  filtering  technique  from  geophysical  data  assimilation  literature,  and  the  Ensemble  Kalman 
Filter  (EnKF).  The  optimized  HHPF  is  found  to  outperform  the  regular  HHPF  using  the  same  particles 
sample  size,  since  the  best  possible  proposal  density  is  used  in  SIS.  For  the  same  particles  sample  sizes, 
the  optimized  HHPF  does  not  perform  as  well  as  the  weights  variance-reduction  particle  filter  and  the 
EnKF  in  terms  of  estimation  accuracy.  By  increasing  the  the  particles  sample  size,  the  optimized  HHPF 
can  achieve  the  same  level  of  accuracy  as  the  other  two  filters.  The  advantage  of  the  optimized  HHPF  is 
in  terms  of  computation  time,  even  with  larger  sample  size  than  the  weights  variance-reduction  particle 
filter  and  EnKF,  due  to  the  incorporation  of  a  homogenization  scheme  that  overcomes  the  need  to  fully 
simulate  the  fast  dynamics.  A  typical  experiment  of  the  optimized  HHPF  using  100  particles,  the  weights 
variance-reduction  particle  filter  and  the  EnKF,  both  using  20  particles,  all  filters  displaying  the  same 
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level  of  accuracy,  has  runtimes  of  134,  1757,  and  540  seconds  for  each  respective  filter.  Some  comparison 
of  the  filter  performances  are  shown  in  the  proceeding  figures. 

Results  of  the  numerical  experiments  has  been  published  [2]. 


Figure  1.  Sparse  observations,  optimized  (Ns  =  100)  and  direct  (Ns  =  400)  HHPFs. 
The  upper  plot  shows  the  estimates  of  an  observed  state,  the  lower  plot  shows  that  of 
an  unobserved  state.  The  unobserved  state  was  estimated  well  by  the  optimized  HHPF 
with  Ns  =  100.  Even  with  Ns  =  400,  The  direct  HHPF  captures  the  fluctuations  in  the 
truth  but  did  not  follow  the  trajectory  well. 


Figure  2.  Sparse  observations,  estimation  error  of  the  observed  states  for  the  optimized 
(Ns  =  100)  and  direct  ( Ns  =  400)  HHPFs  compared  with  observation  error.  The 
optimized  HHPF  estimates  were  as  good  as  observations  (but  the  optimized  HHPF  also 
provided  estimates  of  the  unobserved  states  so  more  information  is  gained  by  using  the 
filter  instead  of  just  observation). 

The  next  step  was  to  consider  the  setting  where  observations  are  sparse  in  time,  i.e.  observations 
are  collected  at  regular  intervals  that  are  large  compared  to  the  numerical  integration  timestep.  In  this 
setting,  for  chaotic  systems,  error  in  between  observation  times  can  grow  significantly  if  the  observation 
interval  is  large  compared  to  the  system’s  error  doubling  time.  Using  the  optimal  control  approach  of 
the  previous  part,  the  problem  can  be  treated  as  a  continuous-time  optimal  control  problem  with  time 
horizon  equal  to  the  length  of  the  observation  interval.  The  cost  function  consists  of  a  running  cost 
that  is  a  quadratic  function  of  the  control  energy,  and  a  terminal  cost  that  is  a  quadratic  function  of 
particle  distance  from  the  location  prescribed  by  the  observation  at  the  end  of  the  observation  interval. 
The  minimum  cost  satisfies  the  Hamilton- Jacobi-Bellman  (HJB)  equation.  The  optimal  control  as  a 
function  of  the  derivative  of  the  minimum  cost,  hence  computation  of  the  optimal  control  requires  solving 
the  HJB  equation.  Using  a  log  transformation,  the  HJB  equation  was  converted  into  a  linear  partial 
differential  equation,  which  solution  can  be  obtained  using  the  Feynman-Kac  formula.  The  derivative 
of  the  Feynman-Kac  representation  involves  the  derivative  of  a  stochastic  path,  which  is  obtained  using 
Malliavin  calculus. 
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Figure  3.  Sparse  observations,  EnKF,  weights  variance-reduction  PF  and  optimized 
HHPF  comparison  at  fixed  Ns  =  100.  The  optimized  HHPF  performed  as  well  as  the 
PF,  but  in  shorter  time  than  both  the  EnKF  and  PF. 
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Figure  4.  Sparse  observations,  comparison  of  estimation  errors  of  the  EnKF,  weights 
variance-reduction  PF  and  optimized  HHPF.  Estimation  error  of  the  optimized  HHPF 
and  the  PF  are  of  the  same  magnitude  at  Ns  =  100.  The  optimized  HHPF  still  required 
less  computation  time  than  the  EnKF  and  PF  even  with  the  EnKF  and  PF  being 
implemented  using  minimum  Ns  possible. 


Augmenting  the  control  modifies  the  signal  dynamics  for  the  particles.  By  the  principle  of  SIS,  this  has 
to  be  accounted  for  by  weighting  the  particles  appropriately  relative  to  the  original  signal  dynamics.  The 
corresponding  weights  are  determined  by  an  application  of  Girsanov’s  Theorem  to  relate  the  probability 
measures  of  the  original  and  controlled  signals.  The  relation  corresponds  to  shifting  the  mean  of  the 
Brownian  motion  in  the  original  signal  by  the  same  amount  as  the  control.  This  weighting  procedure  is 
appended  to  the  particle  weights  based  on  observation  in  the  particle  filtering  algorithm. 

In  practice,  evaluation  of  the  optimal  control  based  on  the  Feynman-Kac  representation  and  Malliavin 
derivative  can  become  computationally  overwhelming  for  nonlinear  signals.  However,  the  optimal  control 
solution  can  be  obtained  explicitly  for  linear  systems.  We  implemented  the  linear  control  strategy  as  a 
suboptimal  control  solution  on  particle  filtering  for  the  3-dimensional  Lorenz  ’63  system.  The  filtering 
results  are  shown  in  Figures  5  and  6.  Implementation  of  the  suboptimal  control  was  sufficient  to  ensure 
consistent  tracking  of  the  signal,  even  when  the  intervals  when  no  observation  is  available  is  large. 
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Figure  5.  Regular  particle  filter;  sparse  observations  collected  at  intervals  equal  to 
error  doubling  time 
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Figure  6.  Particle  filter  with  suboptimal  particle  control;  sparse  observations  collected 
at  intervals  equal  to  error  doubling  time 

2.  Education  and  Academic  Outreach 

Results  from  the  development  of  the  reduced-order  particle  filtering  technique  for  high-dimensional, 
multi-timescale  systems  were  presented  for  the  first  time  at  the  IUTAM  Symposium  on  50  Years  of  Chaos: 
Applied  and  Theoretical  at  Kyoto  University,  Kyoto,  Japan.  The  results  from  an  improved  version  of 
the  algorithm  that  was  capable  of  assimilating  data  on  high-dimensional  chaotic  systems  were  presented 
at  the  IUTAM  Symposium  on  Multiscale  Problems  in  Stochastic  Mechanics  at  Karlsruhe  Institute  of 
Technology,  Karlsruhe,  Germany.  Two  Ph.D  students,  Nishanth  Lingala  and  Hoong  Chieli  Yeong  in 
the  Department  of  Aerospace  Engineering  at  UIUC  were  supported  by  this  grant.  Nicolas  Perkowski 
in  the  Department  of  Mathematics  at  HU,  Berlin,  has  been  working  on  this  research  along  with  his 
graduate  studies  and  was  supported  by  the  grant  during  his  visits  to  UIUC.  Two  other  MS  students, 
Yingtian  Jiang  and  Mahmoud  Mamlouk,  who  have  also  been  involved  part  time,  will  graduate  from  the 
Department  of  Aerospace  Engineering  at  UIUC  in  Fall  2013. 
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